A microscopic model for thin film spreading 
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A microscopic, driven lattice gas model is proposed for the dynamics and spatio-temporal fluc- 
tuations of the precursor film observed in spreading experiments. Matter is transported both by 
holes and particles, and the distribution of each can be described by driven diffusion with a moving 
boundary. This picture leads to a stochastic partial differential equation for the shape of the bound- 
ary, which agrees with the simulations of the lattice gas. Preliminary results for flow in a thermal 
gradient are discussed. 
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Spreading of involatile liquid drops on surfaces and 
fibers plays a significant role in many technologies, the ef- 
ficient deployment of which requires detailed understand- 
ing of the underlying process . Behavior at a macro- 
scopic scale is accurately described by hydrodynamics , 
but experiments performed over the last 15 years or so 
have produced a surprise: in the case of complete wetting, 
the spreading drop, when examined on an atomic length 
scale in the direction normal to the substrate, is found 
to be preceded by a precursor film about one molecule 
thick; this can be followed laterally out to an extension 
of the order of 10^ molecules diameter. At the horizontal 
resolving power of the technique used (ellipsometry), the 
precursor film is found to be flat and homogeneous. Its 
radius advances with time as ^/t. Typical systems show- 
ing this class of behavior which can be investigated by el- 
lipsometry are various silanes spreading over atomically 
flat Si(lll) wafers with highly pure oxydised surfaces. 
With selected silanes, such systems even show dynami- 
cal layering with up to four superposed precursor films 
advancing as y^, the layer directly in contact with the 
substrate being much faster than the others. As Ball ob- 
served over a decade ago |l] , at that time there was hardly 
even a formative theory of such phenomena. Such theory 
must capture the experimentally determined diffusive be- 
havior and, at the same time, explain how an extremely 
viscous, involatile material is transported from the reser- 
voir to the precursor edge. 

A partial solution was first achieved by applying 
Molecular Dynamics to a droplet composed of spheri- 
cal molecules with Lennard- Jones interactions [Q. The 
strength of these is adjusted to achieve sufficient involatil- 
ity. Experimentally, the droplet is placed in contact with 
a substrate composed of much smaller units; this is there- 
fore treated as a continuum for calculating the interac- 
tion of the spreading molecules with the substrate. Such 
a model shows a precursor film with the correct, diffusive 
behavior ||,||. Monte Carlo (MC) lattice gas simulations 
in 3d |^,^| confirm this, and suggest a dual mechanism 
of matter transport involving on the one hand particles 
located in a supernatant layer directly placed above the 



precursor film, and on the other, holes in the precursor 
film itself. At the same time, essentially all the configura- 
tional change occurs at the boundary of the drop. How- 
ever, so far these ideas just constitute a picture without 
quantitative implications. 

In this Letter, we introduce an associated 2d driven 
Ising lattice gas model which can be simulated far more 
effectively. Extensive continuous time MC simulations 
are consistent with the ^/t behavior found in all the ex- 
periments to which we refer |^ . They also allow us to for- 
mulate a probabilistic continuum model of independent 
hole and particle diffusion in the supernatant and precur- 
sor films, which accounts for the matter transfer between 
the reservoir and the precursor edge (PE). The motion 
of the latter in the driven lattice gas is then equivalent 
to a moving boundary in the continuum model. Starting 
from the constitutive equations of the moving boundary 
problem, we derive the effective dynamics of the PE as 
a stochastic partial differential equation. We thus ob- 
tain analytical predictions which are consistent both with 
experiments and with the original 2d Ising model, and 
which can be extended to the asymptotic regime thus far 
unexplored both theoretically and experimentally. 

We now describe our MC lattice gas simulations for 
this problem. At each cell on a square lattice r = (x, y, z) 
with unit lattice spacing, we define the ocuppation num- 
ber n{r,t) to be or 1 depending whether the cell is 
occupied or not. Moreover, we restrict z to the values 
z ~ I (precursor film) and z — 2 (supernatant layer). 
The interaction energy for a given configuration is |q] 



H = -J ^ n{r,t)n{s,t) 

\r-s\ = l 



A J2 ^nir,t), (1) 



where the sums are over r, s on the lattice. The last 
term is the van der Waals interaction with the con- 
tinuum substrate characterized by a Hamaker constant 
A > 0. We are interested in J/ksT large enough to 
achieve a high degree of involatility and A/ksT large 
enough to be in the complete wetting regime Q. Con- 
tinuous time MC Kawasaki dynamics are used |^ with 
initial conditions: (i) n(r,0) — 1 for y = 1, but with 



1 



-s(t) 



X 




y 

FIG. 1. Top view of a typical snapshot of our lattice gas 
model for t — 2 y. lO' MC units. Occupied cells in z — 1 
(precursor film) are in gray, while occupied cells in z = 2 (su- 
pernatant film) appear in black. Non-colored cells are empty. 
Parameters used are A = IQ, J — 1, ksT = 1/3, — 128. 

any x and for both z = 1 and z = 2; (ii) n(r,0) — 
for y > 1. If, as a result of an allowed MC move at 
any later stage, 1, z, t) — 0, then the implied va- 
cancy is filled immediately from the putative reservoir 
droplet, n(x^l^ z,i) = 1 being maintained in this sense 
for all t > 0. Thus, the detailed evolution of the reser- 
voir droplet is ignored. In reality, it would equilibrate 
relatively rapidly in profile by motion of particles near 
its surface and deflate by emission of supernatant par- 
ticles and absorption of precursor holes The main 
advantage of the boundary condition we use, and a very 
significant one, is in computability. We believe our model 
to be new — it is two superposed Ising lattice gases driven 
from one edge. 

A typical snapshot of the simulations is shown in Fig. 
|l|. There is a compact film advancing in both the z — 1 
and z = 2 layers, but the former (shown black) advances 
much faster that the latter (shown gray). It is easy to 
understand why this should be so: if a hole in the precur- 
sor film diffuses to lie beneath the compact film in z = 2, 
at an appropriately selected update it will be filled by 
creating a hole in film above, which diffuses and may be 
trapped at the edge y = 1, or at the supernatant film 
boundary, causing this film to shrink. The precursor ad- 
vances because of emission of holes into itself, but also 
from the arrival of supernatant particles at its edge; such 
particles drop into the z — 1 plane with concomitant 
Hamaker stabilization [ p^ . 

To study the dynamics of the PE, definitions of spin- 
percolative type are needed. Thus, the precursor film 
is the particle cluster P(t) in z = 1 which is con- 
nected along nearest neighbor bonds to the boundary 
line y = 1. The PE is at y = h{x,t) for x — 1, . . . , L^, 
where h{x,t) is given by the maximum value of y among 
those cells {x,y) that belong to V{t). In Fig. || we 
present results for the dynamics of the PE: the veloc- 
ity V{t) = ds/dt [panel (a)], where s{t) is the the mean 
PE displacement s{t) = {^1 Lx)Y^^=ih{.x,t), and the 
roughness u;2(ij,,t) = (C^(a;,t)) [panel (b)]. 
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FIG. 2. Log-log plot of the instantaneous velocity [panel 
(a)] and roughness [panel (b)] of the PE. Solid lines are aver- 
ages over 200 Monte Carlo realizations with ^4 = 10, J = 1, 
fcflT = 1/3, and Lx ~ 256. All units are arbitrary. Dashed 
line in panel (a) depicts the t~^^^ behavior, while those in 
panel (b) correspond to the t^'"^ and t^^^ scaling behaviors. 

where ({x, t) = h{x, t) — s{t) and (• • •) stands for average 
over different MC runs. After a short initial transient 
i = To ~ 10'^ MC units (independent of L^), we recover 
the law s{t) = 7t^/^, with 7 a time independent con- 
stant. This is in clear qualitative agreement with the ex- 
perimental results, and the confirmation of the universal 
y/t law is considerably more precise than in the 3-d simu- 
lations [^ . On the other hand, for times longer than the 
initial transient t > tq — before which the roughness os- 
cillates due to the discreteness of the model and the fiat 
initial condition on the PE — , the PE displays kinetic 
roughening as reflected in the behavior w{L.j.,t) ^ 
[|ll|. The growth exponent /3 seems to cross over from 
the value f3 = 1/6 at intermediate times to (3 = 1/8 
at long times. Due to the restricted range of computa- 
tionally accessible values of L^, the roughness exponent 
a in the stationary state relation w{Lx, 00) ~ L" |Tl| ] 
is more accurately extracted from the behavior of the 
structure factor S{k,t) = {(k{t)C-k{t)) , where Ck{t) is 
the Fourier transform of (^(x,t). For a rough interface, 
S{k,t) is expected [0 to behave as S{k,t) ~ A:-(2"+i) 
for long enough times. In Fig. ^(a) we show the structure 
factor of the PE measured for various times. As is clear 
from the figure, S{k,t) tends to a scale invariant behav- 
ior at long times, characterized by a roughness exponent 
a = 1/2. 

Inspection of typical runs shows a rather sparse distri- 
bution of particles [concentration Cp{r,t)] and holes [con- 
centration Ch{r,t)], which attenuates very rapidly away 
from their respective sources [|o). Since the concentra- 
tions are typically low, the dynamics suggest to introduce 
an independent and non-interacting particle-hole diffu- 
sion model, with identical equations for particles and 
holes. We will treat the r coordinate as a continuous 
variable. Thus, on 7^(i), we have 
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FIG. 3. (a) Structure factor of the PE for t = 1000, 5248, 
29512, 165960, and lO'' MC units (bottom to top). Each 
soHd Hue is an average over 200 MC realizations with A = 10, 
J = 1, ksT = 1/3, and Lx = 512. Units are arbitrary. 
Dashed line is the exact asymptotic structure factor for the 
(discretized) Eq. ([ic|), behaving as S{k, oo) ~ fc"^ for small 
k. (b) Concentration of holes near the PE as a function of 7. 
Dashed line is the analytical prediction (^) and the circles are 
results of MC simulations from fcsT — 0.27 up to fcsT — 0.5, 
left to right. Other parameters are as in panel (a). 



dtCh = DAch - V • g. 



(2) 



where q{r,t) is a noise term representing fluctuations in 
the hole diffusion current within V{t) ||l2|. The interface 
velocity along the outward normal n is given by 



(3) 



The first term on the right in (^) comes from diffusion of 
holes away from the PE. The prcfactor two reflects the 
dual role of particles and holes and that both have the 
same diffusion constant. The second term represents sur- 
face diffusion due to the tendency for particles/holes to 
move along the interface from regions of negative to pos- 
itive curvature, the constant B being proportional to the 
interface concentration of mobile species. Here n is the 
mean curvature. Finally, p{r,t) is a noise term describ- 
ing fluctuations in the interface diffusion current along 
the precursor edge |l^. The concentration obeys also 
the boundary condition Ch(x,0,t) = and the Gibbs- 
Thompson condition Jl4| 



Chir,t)\ 



y=h{x,t) 



Co 



(4) 



where cq is the average concentration of holes at the pre- 
cursor edge and F — caa/kBT, where a is the surface 
tension, which to leading order is taken to be isotropic. 

This moving boundary problem presents formidable 
obstacles to direct attack. Neglecting the noise terms 
q and p, assuming a flat interface, and averaging over x, 
a classical Stefan problem ||l^ is recovered for the mean 
hole concentration c{y,t). On V{t), c{y,t) satisfies the 



diffusion equation dtc = DdyC, but now with far simpler 
boundary conditions c{y,t) \y=g(t) ~ '-O' ^(0, i) — and 

ds 



dt 



(5) 



This problem admits a similarity solution in which 
s{t) = -ft^/'^, and 



7 r^*' 



where 7 is the unique solution of 
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Equation then gives a relation between 7 and cq. 
Now, both 7 and cq can be read off from the simula- 
tions, with different choices of temperature giving differ- 
ent (co,7) points which are compared in Fig. ||(b) with 
the predictions of (|^). The way the simulations of the 
lattice gas are set up implies that, for comparison with 
the diffusive system, we should have D = 1/3. The ex- 
cellent agreement in Fig. ||(b) furnishes very significant 
support for our continuum model. 

The spatio-temporal fluctuations of the PE, expressed 
by w{Lx,t) and S{q,t) in Figs. ||(b) and ||(a), respec- 
tively, require full implementation of the boundary con- 
ditions (||) and a formidable mathematical problem. 
Progress can be made when 7 is small [see Fig. ^(b)] 
by (i) setting dtCh = for y < h{x,t), the quasista- 
tionary approximation p^ ; (ii) perturbing around the 
spatiafly-averaged solution given by (^) and : we 

expand Eq. (|^) to quadratic order in the fluctuation vari- 
able t) = h{x, t) — s{t) describing the precursor edge 
— taking the noise terms to be of the same order as (^{x, t) 
||l2|— , to obtain 



dt 



D\k\ ( rfc2 
v{t) 



Co 

7ii/2 



coih{\k\-ft^''^) + Bk 



(8) 



where Tk [f{x)\ is the Fourier transform of f{x), and r/^ (t) 
is a white noise with correlations [|19l 



{Tlk{t)vk'{t'))=2co5k,-k'5{t-t') 

|fc|27ii/2 



D 



sinh2(|/c|7ii/2) 



B , 



(9) 



In Eq. (g) the linear terms stabilize PE fluctuations, while 
the nonlinear term is of the Kardar-Parisi-Zhang (KPZ) 
type ||2^, with a coefficient V{t) = 7//;^/^ that decays in 
time. In the asymptotic regime i — s- cx), and neglecting 
terms proportional to co, F ~ 0(7^), Eqs. (||) and (||) 
reduce to the conserved linear equation 
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^ = -Bk%+fjk, (10) 

with noise correlations {'nk{t)r)k' {t')) = 2{coB /T)P6k,-k' 
6{t — t'). It is easy to solve this equation exactly and 
obtain (3 = 1/8 and a = 1/2, which provide accurate fits 
at long times for the behavior of the roug hness [Fig. |(b)] 
and the structure factor [Fig. ||(a)], respectively. More- 
over, the behavior of S{k,t) at long but finite times [e.g., 
t = 10^ MC units in Fig. §(a)] is also consistent with 
that predicted by Eq. (|l^). For sufficiently long times, 
the rate of arrival of particles becomes essentially zero. 
In this case, emission of particles to the exterior of the PE 
and holes to the interior are in balance, so the mean inter- 
face is stationary. In such a situation, up to order 0(7^), 
the only fluctuations on the PE are due to the conserved 
surface diffusion current, which is precisely described by 
Eq. (p^). However, for intermediate times {e.g., t < 10^ 
MC units), the behavior of the fluctuations deviates from 
that predicted by Eq. ®, see Figs. |(b) and |(a). A 
plausible explanation is that the KPZ term in Eq. (||) 
is still non-negligible, in which case it should dominate 
the large scale properties of the PE. Using results for 
the KPZ equation | pl] , ^ , the roughness should then be 
given in this time regime hy w ^ [^(i) ^Y^'^ ~ t^^^ , again 
in rather good agreement with the numerical results for 
the lattice gas, see Fig. ^(b). 

In summary, using simulations of our new driven lattice 
gas model we have (i) recaptured the ^/t universal mean 
precursor displacement law; (ii) shown that the spatio- 
temporal fluctuations of the PE demonstrate dynamical 
scaling; (in) shown that the mechanisms of matter trans- 
port from the reservoir droplet to the PE involve low den- 
sity of particles and holes, described by diffusion equa- 
tions with moving boundaries; (iv) derived a stochastic 
partial differential equation for the PE which brings in 
the fluctuation regimes found previously in a natural way, 
and which tends asymptotically to a Langevin descrip- 
tion for a free interface with conserved dynamics as in 
(|lO|), possibly observable experimentally, even though it 
presents a KPZ behavior at intermediate times. In pre- 
liminary work, we have investigated our lattice gas model 
for a single layer (z = 1) in a thermal gradient, finding a 
mean spatial displacement growing linearly in time, gen- 
erated by biased hole diffusion. This has relevance for 
microfluidics 
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